Modeling relative establishment probability
library(here)## here() starts at /Users/julienvollering/Desktop/PortableOffice/conifer-plantation-spreadlibrary(tidyverse)## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang## ── Attaching packages ────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──## ✔ ggplot2 3.1.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0## ── Conflicts ───────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(glmmTMB)
library(DHARMa)## Registered S3 method overwritten by 'DHARMa':
## method from
## refit.glmmTMB glmmTMBdatlist <- readRDS(here("output","data.rds"))
locality.preds <- c('age','bio01','bio19','seeds.ExP','seeds.WALD','relelev')
contrast <- "T4"
excluded <- "plantation forest"P.sitchensis-lutz
Ps_dat <- Filter(function(x)!all(is.na(x)), datlist$`P.sitchensis-lutz`)Ps_arith <- pivot_longer(Ps_dat, 10:ncol(Ps_dat), names_to = "nin", values_to = "are") %>%
group_by(nin) %>%
summarize(daa = sum(are)/10, # units from are to decare
n = sum(wildlings*are)) %>%
mutate(density = n/daa) %>%
arrange(desc(density), desc(n), daa)
scaler <- filter(Ps_arith, nin == contrast) %>%
pull(density)
Ps_arith <- mutate(Ps_arith, fit = density/scaler)| nin | daa | n | density | fit |
|---|---|---|---|---|
| T17 | 2.0 | 42.8 | 21.1 | 5.4 |
| T2 | 160.2 | 961.2 | 6.0 | 1.5 |
| T34 | 1071.0 | 5664.4 | 5.3 | 1.3 |
| T4 | 1362.9 | 5360.0 | 3.9 | 1.0 |
| fen | 758.3 | 2886.0 | 3.8 | 1.0 |
| disturbed | 409.5 | 1550.0 | 3.8 | 1.0 |
| plantation forest | 88.9 | 227.5 | 2.6 | 0.7 |
| T32 | 676.5 | 1236.0 | 1.8 | 0.5 |
| T16 | 14.8 | 25.3 | 1.7 | 0.4 |
| swamp | 160.5 | 251.0 | 1.6 | 0.4 |
| bog | 186.8 | 252.8 | 1.4 | 0.3 |
| T6 | 34.9 | 39.4 | 1.1 | 0.3 |
| T1 | 141.5 | 158.0 | 1.1 | 0.3 |
| T27 | 68.1 | 75.3 | 1.1 | 0.3 |
| T31 | 317.4 | 342.3 | 1.1 | 0.3 |
| T21 | 6.3 | 6.4 | 1.0 | 0.3 |
| T29 | 6.7 | 4.8 | 0.7 | 0.2 |
| cultivated | 1830.1 | 1015.8 | 0.6 | 0.1 |
| T24 | 6.9 | 3.1 | 0.5 | 0.1 |
| T12 | 11.7 | 4.6 | 0.4 | 0.1 |
| T13 | 11.7 | 2.0 | 0.2 | 0.0 |
| T18 | 2.0 | 0.0 | 0.0 | 0.0 |
| T3 | 2.3 | 0.0 | 0.0 | 0.0 |
| T33 | 5.9 | 0.0 | 0.0 | 0.0 |
| spring | 9.8 | 0.0 | 0.0 | 0.0 |
perfectlyseparated <- c('T18', 'T3', 'T33', 'spring')
removed <- c(perfectlyseparated, excluded)
Ps_dat <- Ps_dat %>%
filter_at(removed, all_vars(. < 0.5)) %>%
filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>%
select(-c('species'))Ps_means <- Ps_dat[, locality.preds] %>% map_dbl(mean)
Ps_sds <- Ps_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
Ps_dat <- modify_at(Ps_dat, .at = i, ~ (. - Ps_means[i])/Ps_sds[i])
}
summary(Ps_dat)## locality wildlings age bio01
## Length:73309 Min. : 0.0000 Min. :-1.5919 Min. :-2.3687
## Class :character 1st Qu.: 0.0000 1st Qu.:-0.9296 1st Qu.:-0.5921
## Mode :character Median : 0.0000 Median : 0.1006 Median :-0.3229
## Mean : 0.2788 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.0000 3rd Qu.: 0.5422 3rd Qu.: 0.9082
## Max. :206.0000 Max. : 2.2347 Max. : 1.5957
## bio19 seeds.ExP seeds.WALD relelev
## Min. :-1.8113 Min. :-0.2245 Min. :-0.49584 Min. :-3.3307
## 1st Qu.:-0.5378 1st Qu.:-0.2245 1st Qu.:-0.47446 1st Qu.:-0.4699
## Median :-0.2429 Median :-0.2245 Median :-0.40001 Median : 0.0833
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.9368 3rd Qu.:-0.2241 3rd Qu.:-0.05296 3rd Qu.: 0.4406
## Max. : 2.5052 Max. :13.1032 Max. :14.58650 Max. : 7.9503
## plantation forest cultivated disturbed T32
## Min. :0.0000000 Min. :0.0000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.0000000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.0000000 Median :0.0000 Median :0.00000 Median :0.00000
## Mean :0.0006327 Mean :0.2495 Mean :0.05579 Mean :0.09212
## 3rd Qu.:0.0000000 3rd Qu.:0.5250 3rd Qu.:0.00000 3rd Qu.:0.00000
## Max. :0.4910000 Max. :1.0000 Max. :1.00000 Max. :1.00000
## fen T34 T13 T4
## Min. :0.0000 Min. :0.0000 Min. :0.000000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000000 1st Qu.:0.0000
## Median :0.0000 Median :0.0000 Median :0.000000 Median :0.0000
## Mean :0.1033 Mean :0.1458 Mean :0.001598 Mean :0.1858
## 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.000000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :1.000000 Max. :1.0000
## T2 T1 swamp T18
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.00e+00
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.00e+00
## Median :0.00000 Median :0.00000 Median :0.0000 Median :0.00e+00
## Mean :0.02185 Mean :0.01929 Mean :0.0219 Mean :7.64e-05
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.00e+00
## Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :4.67e-01
## T27 T17 spring
## Min. :0.000000 Min. :0.0000000 Min. :0.000000
## 1st Qu.:0.000000 1st Qu.:0.0000000 1st Qu.:0.000000
## Median :0.000000 Median :0.0000000 Median :0.000000
## Mean :0.009282 Mean :0.0002765 Mean :0.001233
## 3rd Qu.:0.000000 3rd Qu.:0.0000000 3rd Qu.:0.000000
## Max. :0.999000 Max. :1.0000000 Max. :0.499000
## bog T6 T31 T16
## Min. :0.00000 Min. :0.000000 Min. :0.00000 Min. :0.000000
## 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.000000
## Median :0.00000 Median :0.000000 Median :0.00000 Median :0.000000
## Mean :0.02548 Mean :0.004763 Mean :0.04324 Mean :0.002025
## 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.:0.00000 3rd Qu.:0.000000
## Max. :1.00000 Max. :1.000000 Max. :1.00000 Max. :1.000000
## T29 T12 T24
## Min. :0.000000 Min. :0.000000 Min. :0.000000
## 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.000000
## Median :0.000000 Median :0.000000 Median :0.000000
## Mean :0.000916 Mean :0.001592 Mean :0.000945
## 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.:0.000000
## Max. :0.922000 Max. :1.000000 Max. :0.500000
## T3 T21 T33
## Min. :0.00e+00 Min. :0.0000000 Min. :0.000000
## 1st Qu.:0.00e+00 1st Qu.:0.0000000 1st Qu.:0.000000
## Median :0.00e+00 Median :0.0000000 Median :0.000000
## Mean :4.52e-05 Mean :0.0008594 Mean :0.000216
## 3rd Qu.:0.00e+00 3rd Qu.:0.0000000 3rd Qu.:0.000000
## Max. :4.27e-01 Max. :1.0000000 Max. :0.498000paste(names(Ps_dat)[!(names(Ps_dat) %in% c(contrast, removed))], collapse = " + ")## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21"Compare WALD vs ExP
Ps_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 +
swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 +
T21 + (1 | locality)")
Ps_m1 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ 0)
summary(Ps_m1)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +
## T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |
## locality)
## Data: Ps_dat
##
## AIC BIC logLik deviance df.resid
## 44876.3 45133.9 -22410.1 44820.3 73281
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 3.745 1.935
## Number of obs: 73309, groups: locality, 36
##
## Overdispersion parameter for nbinom2 family (): 0.0788
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3588408 0.3312308 -7.12 1.07e-12 ***
## seeds.WALD 0.9384915 0.0266383 35.23 < 2e-16 ***
## age 0.2583634 0.3360831 0.77 0.442043
## bio01 0.5949349 0.4141349 1.44 0.150839
## bio19 -0.5308042 0.4021168 -1.32 0.186827
## relelev -0.2617019 0.0360642 -7.26 3.97e-13 ***
## cultivated -2.6215100 0.1013181 -25.87 < 2e-16 ***
## disturbed -0.1931504 0.1223341 -1.58 0.114364
## T32 -0.9567007 0.1188110 -8.05 8.13e-16 ***
## fen -0.4023827 0.1093915 -3.68 0.000235 ***
## T34 -0.1604257 0.1097421 -1.46 0.143784
## T13 0.3284157 0.8961519 0.37 0.714012
## T2 0.0008741 0.1998395 0.00 0.996510
## T1 -2.0711687 0.2804754 -7.38 1.53e-13 ***
## swamp -2.3173703 0.2126976 -10.90 < 2e-16 ***
## T27 1.1609938 0.4601492 2.52 0.011633 *
## T17 2.8532089 1.0404665 2.74 0.006102 **
## bog -1.5504993 0.1971566 -7.86 3.71e-15 ***
## T6 -1.4496908 0.3818381 -3.80 0.000147 ***
## T31 1.0036357 0.2054121 4.89 1.03e-06 ***
## T16 0.3939334 0.4827710 0.82 0.414509
## T29 -1.4845290 1.1672067 -1.27 0.203421
## T12 -1.8898514 1.4458260 -1.31 0.191176
## T24 -3.1043266 2.1623452 -1.44 0.151108
## T21 -3.4083792 0.7636041 -4.46 8.06e-06 ***
## seeds.WALD:age 0.1599230 0.0259332 6.17 6.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ps_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T34 + T13 + T2 + T1 +
swamp + T27 + T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 +
T21 + (1 | locality)")
Ps_m2 <- glmmTMB(Ps_f2, Ps_dat, family = "nbinom2", ziformula = ~ 0)
summary(Ps_m2)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +
## T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |
## locality)
## Data: Ps_dat
##
## AIC BIC logLik deviance df.resid
## 46062.6 46320.2 -23003.3 46006.6 73281
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 3.268 1.808
## Number of obs: 73309, groups: locality, 36
##
## Overdispersion parameter for nbinom2 family (): 0.0639
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.00098 0.31012 -6.452 1.10e-10 ***
## seeds.ExP 0.52240 0.02614 19.985 < 2e-16 ***
## age 0.40449 0.31396 1.288 0.197629
## bio01 0.50619 0.38742 1.307 0.191359
## bio19 -0.59032 0.37635 -1.569 0.116761
## relelev -0.31592 0.03810 -8.292 < 2e-16 ***
## cultivated -2.90729 0.10156 -28.625 < 2e-16 ***
## disturbed -0.12715 0.12970 -0.980 0.326947
## T32 -0.98601 0.11678 -8.444 < 2e-16 ***
## fen -0.61490 0.11376 -5.405 6.47e-08 ***
## T34 -0.15135 0.11279 -1.342 0.179644
## T13 0.14798 0.92685 0.160 0.873153
## T2 0.29751 0.21278 1.398 0.162057
## T1 -2.77607 0.28435 -9.763 < 2e-16 ***
## swamp -2.62667 0.22223 -11.819 < 2e-16 ***
## T27 0.74253 0.50100 1.482 0.138316
## T17 3.22023 1.19005 2.706 0.006811 **
## bog -1.34171 0.19957 -6.723 1.78e-11 ***
## T6 -1.30451 0.40234 -3.242 0.001186 **
## T31 1.05021 0.20777 5.055 4.31e-07 ***
## T16 0.43948 0.51149 0.859 0.390216
## T29 -0.81335 1.20018 -0.678 0.497965
## T12 -2.57557 1.44967 -1.777 0.075624 .
## T24 -0.67897 1.88438 -0.360 0.718614
## T21 -2.71096 0.76065 -3.564 0.000365 ***
## seeds.ExP:age -0.03212 0.02600 -1.235 0.216686
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Ps_m1 0.0 28
## Ps_m2 1186.3 28Continue with WALD.
Compare zero inflation effects
Ps_m3 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ 1)
summary(Ps_m3)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +
## T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |
## locality)
## Zero inflation: ~1
## Data: Ps_dat
##
## AIC BIC logLik deviance df.resid
## 44878.3 45145.1 -22410.1 44820.3 73280
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 3.745 1.935
## Number of obs: 73309, groups: locality, 36
##
## Overdispersion parameter for nbinom2 family (): 0.0788
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3588826 0.3312219 -7.12 1.07e-12 ***
## seeds.WALD 0.9385074 0.0266387 35.23 < 2e-16 ***
## age 0.2582882 0.3360732 0.77 0.442162
## bio01 0.5948216 0.4141232 1.44 0.150906
## bio19 -0.5308007 0.4021053 -1.32 0.186817
## relelev -0.2616964 0.0360642 -7.26 3.98e-13 ***
## cultivated -2.6215078 0.1013185 -25.87 < 2e-16 ***
## disturbed -0.1932046 0.1223341 -1.58 0.114263
## T32 -0.9566731 0.1188114 -8.05 8.14e-16 ***
## fen -0.4023916 0.1093918 -3.68 0.000235 ***
## T34 -0.1604635 0.1097424 -1.46 0.143691
## T13 0.3284436 0.8961246 0.37 0.713980
## T2 0.0008862 0.1998404 0.00 0.996462
## T1 -2.0711986 0.2804756 -7.38 1.53e-13 ***
## swamp -2.3173310 0.2126978 -10.89 < 2e-16 ***
## T27 1.1607030 0.4601308 2.52 0.011651 *
## T17 2.8532586 1.0404906 2.74 0.006102 **
## bog -1.5504769 0.1971573 -7.86 3.72e-15 ***
## T6 -1.4496469 0.3818395 -3.80 0.000147 ***
## T31 1.0035677 0.2054108 4.89 1.03e-06 ***
## T16 0.3940391 0.4827768 0.82 0.414390
## T29 -1.4845403 1.1672111 -1.27 0.203419
## T12 -1.8896773 1.4457986 -1.31 0.191208
## T24 -3.1045989 2.1623385 -1.44 0.151071
## T21 -3.4085395 0.7636025 -4.46 8.05e-06 ***
## seeds.WALD:age 0.1599326 0.0259335 6.17 6.96e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.33 550.26 -0.031 0.975Ps_m4 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ age)
summary(Ps_m4)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +
## T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |
## locality)
## Zero inflation: ~age
## Data: Ps_dat
##
## AIC BIC logLik deviance df.resid
## 44355.8 44631.9 -22147.9 44295.8 73279
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 3.744 1.935
## Number of obs: 73309, groups: locality, 36
##
## Overdispersion parameter for nbinom2 family (): 0.16
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.51775 0.33423 -4.54 5.60e-06 ***
## seeds.WALD 0.90437 0.02792 32.40 < 2e-16 ***
## age -0.25650 0.33716 -0.76 0.446792
## bio01 0.61522 0.41432 1.48 0.137576
## bio19 -0.53728 0.40203 -1.34 0.181411
## relelev -0.26289 0.03411 -7.71 1.28e-14 ***
## cultivated -2.58291 0.10399 -24.84 < 2e-16 ***
## disturbed -0.43691 0.12949 -3.37 0.000741 ***
## T32 -1.09632 0.11718 -9.36 < 2e-16 ***
## fen -0.39970 0.11135 -3.59 0.000331 ***
## T34 -0.14437 0.10765 -1.34 0.179893
## T13 0.32869 0.89216 0.37 0.712565
## T2 -0.87796 0.24094 -3.64 0.000269 ***
## T1 -2.34639 0.31287 -7.50 6.40e-14 ***
## swamp -2.38736 0.25218 -9.47 < 2e-16 ***
## T27 1.42264 0.48746 2.92 0.003517 **
## T17 2.50901 0.77584 3.23 0.001221 **
## bog -1.70283 0.22673 -7.51 5.90e-14 ***
## T6 -2.03677 0.39124 -5.21 1.93e-07 ***
## T31 1.03591 0.20327 5.10 3.47e-07 ***
## T16 0.40353 0.45464 0.89 0.374768
## T29 -2.16010 1.48667 -1.45 0.146228
## T12 -2.35583 1.83304 -1.29 0.198723
## T24 -3.34399 2.79879 -1.19 0.232165
## T21 -3.28725 0.72425 -4.54 5.66e-06 ***
## seeds.WALD:age 0.11690 0.02579 4.53 5.84e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.14693 0.09179 -1.601 0.109
## age -1.11316 0.05366 -20.747 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ps_m5 <- glmmTMB(Ps_f1, Ps_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Ps_m5)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T13 + T2 + T1 + swamp + T27 +
## T17 + bog + T6 + T31 + T16 + T29 + T12 + T24 + T21 + (1 |
## locality)
## Zero inflation: ~age + (1 | locality)
## Data: Ps_dat
##
## AIC BIC logLik deviance df.resid
## 42266.6 42551.9 -21102.3 42204.6 73278
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 2.034 1.426
## Number of obs: 73309, groups: locality, 36
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 5.133 2.266
## Number of obs: 73309, groups: locality, 36
##
## Overdispersion parameter for nbinom2 family (): 0.274
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.35045 0.26151 -1.34 0.180208
## seeds.WALD 0.83919 0.02646 31.72 < 2e-16 ***
## age 0.04685 0.26979 0.17 0.862133
## bio01 -0.15971 0.40960 -0.39 0.696587
## bio19 -0.86385 0.38061 -2.27 0.023230 *
## relelev -0.30122 0.03738 -8.06 7.72e-16 ***
## cultivated -2.63220 0.10448 -25.19 < 2e-16 ***
## disturbed -0.76252 0.15240 -5.00 5.63e-07 ***
## T32 -0.80398 0.11738 -6.85 7.40e-12 ***
## fen -0.15901 0.11135 -1.43 0.153271
## T34 0.03014 0.10599 0.28 0.776104
## T13 0.47073 0.91836 0.51 0.608249
## T2 1.04641 0.26711 3.92 8.95e-05 ***
## T1 -2.32908 0.24152 -9.64 < 2e-16 ***
## swamp -2.02971 0.29586 -6.86 6.86e-12 ***
## T27 1.56716 0.42690 3.67 0.000242 ***
## T17 2.30316 0.63314 3.64 0.000275 ***
## bog -1.35952 0.27339 -4.97 6.60e-07 ***
## T6 0.72195 0.74465 0.97 0.332290
## T31 0.73905 0.23772 3.11 0.001877 **
## T16 0.51663 0.50709 1.02 0.308294
## T29 -2.27767 1.85349 -1.23 0.219127
## T12 -2.79530 1.87336 -1.49 0.135664
## T24 -1.86526 2.58922 -0.72 0.471283
## T21 -5.10913 1.31469 -3.89 0.000102 ***
## seeds.WALD:age 0.18578 0.02372 7.83 4.78e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2857 0.3924 3.277 0.00105 **
## age -0.8856 0.3880 -2.282 0.02247 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Ps_m5 0.0 31
## Ps_m4 2089.2 30
## Ps_m1 2609.6 28
## Ps_m3 2611.6 29m5 is better than m1, m3, and m4: much better AIC, clear zero-inflation effects
simulationOutput <- simulateResiduals(fittedModel = Ps_m5)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.075717, p-value = 0.104
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0202, p-value = 0.248
## alternative hypothesis: two.sidedDiagnostics look OK.
Interpret relative establishment probability
selmod <- Ps_m5
terms <- attributes(selmod$frame)$terms
newdat <- attributes(terms)$factors %>%
as_tibble() %>%
select(-`seeds.WALD:age`) %>%
mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Ps_means["seeds.WALD"]) / Ps_sds["seeds.WALD"],
age = (Ps_means["age"] - Ps_means["age"]) / Ps_sds["age"],
bio01 = (Ps_means["bio01"] - Ps_means["bio01"]) / Ps_sds["bio01"],
bio19 = (Ps_means["bio19"] - Ps_means["bio19"]) / Ps_sds["bio19"],
relelev = (Ps_means["relelev"] - Ps_means["relelev"]) / Ps_sds["relelev"],
locality = NA_character_) %>%
distinct()
newdat <- newdat %>%
mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>%
replace_na(list(nin = "T4"))Ps_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Ps_link <- mutate(Ps_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>%
select(-se.fit) %>%
mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Ps_link, nin == "T4") %>%
pull(fit)
link.scaled <- mutate_at(Ps_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) +
geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) +
geom_point() +
geom_point(data = filter(Ps_arith, nin %in% link.scaled$nin),
col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()
newdat.locality <- newdat %>%
select(-locality) %>%
expand_grid(locality = c(NA_character_, unique(Ps_dat$locality)))
resp.locality <- newdat.locality %>%
bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>%
pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))
gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) +
geom_point(data = filter(resp.locality.scaled, !is.na(locality)),
mapping = aes(color = locality), cex = 0.5, show.legend = FALSE) +
geom_point(data = filter(resp.locality.scaled, is.na(locality))) +
geom_point(data = filter(Ps_arith, nin %in% resp.locality.scaled$nin),
color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()
Larix
L_dat <- Filter(function(x)!all(is.na(x)), datlist$Larix)L_arith <- pivot_longer(L_dat, 10:ncol(L_dat), names_to = "nin", values_to = "are") %>%
group_by(nin) %>%
summarize(daa = sum(are)/10, # units from are to decare
n = sum(wildlings*are)) %>%
mutate(density = n/daa) %>%
arrange(desc(density), desc(n), daa)
scaler <- filter(L_arith, nin == contrast) %>%
pull(density)
L_arith <- mutate(L_arith, fit = density/scaler)| nin | daa | n | density | fit |
|---|---|---|---|---|
| T17 | 2.0 | 211.8 | 104.5 | 160.6 |
| plantation forest | 169.1 | 159.8 | 0.9 | 1.5 |
| T34 | 471.2 | 308.7 | 0.7 | 1.0 |
| T4 | 332.4 | 216.3 | 0.7 | 1.0 |
| T32 | 227.5 | 115.0 | 0.5 | 0.8 |
| T1 | 40.3 | 6.9 | 0.2 | 0.3 |
| disturbed | 185.8 | 17.5 | 0.1 | 0.1 |
| cultivated | 1243.6 | 91.1 | 0.1 | 0.1 |
| fen | 132.8 | 6.1 | 0.0 | 0.1 |
| T18 | 0.6 | 0.0 | 0.0 | 0.0 |
| T30 | 1.0 | 0.0 | 0.0 | 0.0 |
| T27 | 1.5 | 0.0 | 0.0 | 0.0 |
| T2 | 1.6 | 0.0 | 0.0 | 0.0 |
| swamp | 1.7 | 0.0 | 0.0 | 0.0 |
| T13 | 3.0 | 0.0 | 0.0 | 0.0 |
perfectlyseparated <- c('T18', 'T30', 'T2', 'T27', 'swamp', 'T13')
removed <- c(perfectlyseparated, excluded)
L_dat <- L_dat %>%
filter_at(removed, all_vars(. < 0.5)) %>%
filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>%
select(-c('species'))L_means <- L_dat[, locality.preds] %>% map_dbl(mean)
L_sds <- L_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
L_dat <- modify_at(L_dat, .at = i, ~ (. - L_means[i])/L_sds[i])
}
summary(L_dat)## locality wildlings age
## Length:26514 Min. : 0.00000 Min. :-1.9249
## Class :character 1st Qu.: 0.00000 1st Qu.:-0.2941
## Mode :character Median : 0.00000 Median : 0.3353
## Mean : 0.03719 Mean : 0.0000
## 3rd Qu.: 0.00000 3rd Qu.: 0.5642
## Max. :27.00000 Max. : 1.5655
## bio01 bio19 seeds.ExP seeds.WALD
## Min. :-2.53170 Min. :-1.5815 Min. :-0.1189 Min. :-0.3842
## 1st Qu.: 0.05949 1st Qu.:-0.4381 1st Qu.:-0.1189 1st Qu.:-0.3666
## Median : 0.29676 Median :-0.2804 Median :-0.1189 Median :-0.3109
## Mean : 0.00000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.67912 3rd Qu.: 0.4687 3rd Qu.:-0.1189 3rd Qu.:-0.1421
## Max. : 0.90390 Max. : 2.2034 Max. :21.3893 Max. :19.1019
## relelev plantation forest cultivated disturbed
## Min. :-3.16949 Min. :0.00000 Min. :0.0000 Min. :0.00000
## 1st Qu.:-0.54556 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.00000
## Median :-0.03968 Median :0.00000 Median :0.1405 Median :0.00000
## Mean : 0.00000 Mean :0.00384 Mean :0.4678 Mean :0.06954
## 3rd Qu.: 0.46001 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.00000
## Max. : 5.73614 Max. :0.49900 Max. :1.0000 Max. :1.00000
## T32 fen T34 T13
## Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.0000000
## 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000000
## Median :0.00000 Median :0.00000 Median :0.0000 Median :0.0000000
## Mean :0.08533 Mean :0.04989 Mean :0.1766 Mean :0.0002235
## 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000000
## Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :0.4650000
## T4 T2 T1 swamp
## Min. :0.0000 Min. :0.0000000 Min. :0.00000 Min. :0.00e+00
## 1st Qu.:0.0000 1st Qu.:0.0000000 1st Qu.:0.00000 1st Qu.:0.00e+00
## Median :0.0000 Median :0.0000000 Median :0.00000 Median :0.00e+00
## Mean :0.1247 Mean :0.0002907 Mean :0.01493 Mean :6.74e-05
## 3rd Qu.:0.0000 3rd Qu.:0.0000000 3rd Qu.:0.00000 3rd Qu.:0.00e+00
## Max. :1.0000 Max. :0.4980000 Max. :1.00000 Max. :4.83e-01
## T30 T18 T27 T17
## Min. :0 Min. :0 Min. :0.000000 Min. :0.0000000
## 1st Qu.:0 1st Qu.:0 1st Qu.:0.000000 1st Qu.:0.0000000
## Median :0 Median :0 Median :0.000000 Median :0.0000000
## Mean :0 Mean :0 Mean :0.000159 Mean :0.0007646
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0.000000 3rd Qu.:0.0000000
## Max. :0 Max. :0 Max. :0.419000 Max. :1.0000000paste(names(L_dat)[!(names(L_dat) %in% c(contrast, removed))], collapse = " + ")## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T34 + T1 + T17"Compare WALD vs ExP
L_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T34 + T1 + T17 +
(1 | locality)")
L_m1 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ 0)
summary(L_m1)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Data: L_dat
##
## AIC BIC logLik deviance df.resid
## 3878.1 4009.1 -1923.1 3846.1 26498
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2738 0.5233
## Number of obs: 26514, groups: locality, 12
##
## Overdispersion parameter for nbinom2 family (): 0.0595
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.84395 0.29882 -12.864 < 2e-16 ***
## seeds.WALD 1.03433 0.06119 16.903 < 2e-16 ***
## age -0.44714 0.31379 -1.425 0.154170
## bio01 -1.03506 0.39701 -2.607 0.009130 **
## bio19 -1.93007 0.47899 -4.029 5.59e-05 ***
## relelev 0.02552 0.08020 0.318 0.750277
## cultivated -2.43001 0.32461 -7.486 7.10e-14 ***
## disturbed -2.41780 0.46959 -5.149 2.62e-07 ***
## T32 -1.22154 0.32713 -3.734 0.000188 ***
## fen -3.72529 0.83550 -4.459 8.24e-06 ***
## T34 0.03735 0.35197 0.106 0.915484
## T1 -1.94085 0.97643 -1.988 0.046844 *
## T17 6.40605 1.05478 6.073 1.25e-09 ***
## seeds.WALD:age 0.24491 0.06875 3.562 0.000367 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1L_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T34 + T1 + T17 +
(1 | locality)")
L_m2 <- glmmTMB(L_f2, L_dat, family = "nbinom2", ziformula = ~ 0)
summary(L_m2)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Data: L_dat
##
## AIC BIC logLik deviance df.resid
## 4249.0 4379.9 -2108.5 4217.0 26498
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.7066 0.8406
## Number of obs: 26514, groups: locality, 12
##
## Overdispersion parameter for nbinom2 family (): 0.0331
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.42176 0.35463 -6.829 8.55e-12 ***
## seeds.ExP 0.64949 0.08092 8.026 1.01e-15 ***
## age 0.41651 0.43886 0.949 0.342584
## bio01 0.21247 0.56587 0.375 0.707303
## bio19 -0.56967 0.67811 -0.840 0.400861
## relelev -0.03232 0.09079 -0.356 0.721811
## cultivated -3.63714 0.34925 -10.414 < 2e-16 ***
## disturbed -3.27412 0.46437 -7.051 1.78e-12 ***
## T32 -1.72551 0.34102 -5.060 4.20e-07 ***
## fen -4.50569 0.65557 -6.873 6.29e-12 ***
## T34 -1.36713 0.37758 -3.621 0.000294 ***
## T1 -1.92174 0.82430 -2.331 0.019734 *
## T17 5.04949 1.33503 3.782 0.000155 ***
## seeds.ExP:age -0.07588 0.08216 -0.924 0.355696
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## L_m1 0.0 16
## L_m2 370.8 16Continue with WALD.
Compare zero inflation effects
L_m3 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ 1)
summary(L_m3)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation: ~1
## Data: L_dat
##
## AIC BIC logLik deviance df.resid
## 3880.1 4019.3 -1923.1 3846.1 26497
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2738 0.5233
## Number of obs: 26514, groups: locality, 12
##
## Overdispersion parameter for nbinom2 family (): 0.0595
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.84392 0.29882 -12.864 < 2e-16 ***
## seeds.WALD 1.03432 0.06119 16.903 < 2e-16 ***
## age -0.44713 0.31379 -1.425 0.154181
## bio01 -1.03507 0.39702 -2.607 0.009130 **
## bio19 -1.93010 0.47900 -4.029 5.59e-05 ***
## relelev 0.02552 0.08020 0.318 0.750276
## cultivated -2.43007 0.32461 -7.486 7.10e-14 ***
## disturbed -2.41785 0.46959 -5.149 2.62e-07 ***
## T32 -1.22158 0.32713 -3.734 0.000188 ***
## fen -3.72532 0.83550 -4.459 8.24e-06 ***
## T34 0.03731 0.35197 0.106 0.915583
## T1 -1.94086 0.97642 -1.988 0.046843 *
## T17 6.40607 1.05479 6.073 1.25e-09 ***
## seeds.WALD:age 0.24490 0.06875 3.562 0.000367 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.7 2484.9 -0.007 0.995L_m4 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ age)
summary(L_m4)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation: ~age
## Data: L_dat
##
## AIC BIC logLik deviance df.resid
## 3869.9 4017.3 -1917.0 3833.9 26496
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2046 0.4523
## Number of obs: 26514, groups: locality, 12
##
## Overdispersion parameter for nbinom2 family (): 0.0644
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.52811 0.30007 -11.758 < 2e-16 ***
## seeds.WALD 1.02067 0.06336 16.109 < 2e-16 ***
## age -0.98463 0.32757 -3.006 0.00265 **
## bio01 -1.19452 0.36823 -3.244 0.00118 **
## bio19 -2.06985 0.43459 -4.763 1.91e-06 ***
## relelev 0.06530 0.09487 0.688 0.49125
## cultivated -2.42125 0.31614 -7.659 1.88e-14 ***
## disturbed -2.43422 0.46514 -5.233 1.66e-07 ***
## T32 -1.26467 0.32281 -3.918 8.94e-05 ***
## fen -3.79127 0.82779 -4.580 4.65e-06 ***
## T34 -0.02954 0.34640 -0.085 0.93204
## T1 -1.86342 0.94799 -1.966 0.04934 *
## T17 6.39796 1.01722 6.290 3.18e-10 ***
## seeds.WALD:age 0.25390 0.09144 2.777 0.00549 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.400 2.975 -1.143 0.253
## age -2.458 1.537 -1.599 0.110L_m5 <- glmmTMB(L_f1, L_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(L_m5)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T34 + T1 + T17 + (1 | locality)
## Zero inflation: ~age + (1 | locality)
## Data: L_dat
##
## AIC BIC logLik deviance df.resid
## 3843.2 3998.7 -1902.6 3805.2 26495
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2603 0.5102
## Number of obs: 26514, groups: locality, 12
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 1.411 1.188
## Number of obs: 26514, groups: locality, 12
##
## Overdispersion parameter for nbinom2 family (): 0.138
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.77051 0.36671 -7.555 4.19e-14 ***
## seeds.WALD 1.04741 0.06627 15.804 < 2e-16 ***
## age -1.23721 0.42699 -2.898 0.00376 **
## bio01 -0.96633 0.43191 -2.237 0.02526 *
## bio19 -2.07193 0.52678 -3.933 8.38e-05 ***
## relelev 0.04995 0.09175 0.544 0.58612
## cultivated -2.45366 0.32924 -7.452 9.17e-14 ***
## disturbed -2.39472 0.48467 -4.941 7.78e-07 ***
## T32 -1.37959 0.33335 -4.139 3.49e-05 ***
## fen -3.53032 0.80225 -4.401 1.08e-05 ***
## T34 0.14793 0.36357 0.407 0.68409
## T1 -2.02332 0.87769 -2.305 0.02115 *
## T17 6.05552 0.71886 8.424 < 2e-16 ***
## seeds.WALD:age 0.26766 0.10735 2.493 0.01265 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4100 0.5420 0.756 0.4494
## age -1.0206 0.4953 -2.061 0.0393 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## L_m5 0.0 19
## L_m4 26.7 18
## L_m1 34.9 16
## L_m3 36.9 17m5 is better than m1, m3, and m4: much better AIC, clear zero-inflation effects
simulationOutput <- simulateResiduals(fittedModel = L_m5)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.25346, p-value = 0.144
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99942, p-value = 0.816
## alternative hypothesis: two.sidedDiagnostics look OK.
Interpret relative establishment probability
selmod <- L_m5
terms <- attributes(selmod$frame)$terms
newdat <- attributes(terms)$factors %>%
as_tibble() %>%
select(-`seeds.WALD:age`) %>%
mutate(seeds.WALD = (Ps_means["seeds.WALD"] - L_means["seeds.WALD"]) / L_sds["seeds.WALD"],
age = (Ps_means["age"] - L_means["age"]) / L_sds["age"],
bio01 = (Ps_means["bio01"] - L_means["bio01"]) / L_sds["bio01"],
bio19 = (Ps_means["bio19"] - L_means["bio19"]) / L_sds["bio19"],
relelev = (Ps_means["relelev"] - L_means["relelev"]) / L_sds["relelev"],
locality = NA_character_) %>%
distinct()
newdat <- newdat %>%
mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>%
replace_na(list(nin = "T4"))L_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
L_link <- mutate(L_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>%
select(-se.fit) %>%
mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(L_link, nin == "T4") %>%
pull(fit)
link.scaled <- mutate_at(L_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) +
geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) +
geom_point() +
geom_point(data = filter(L_arith, nin %in% link.scaled$nin),
col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()
newdat.locality <- newdat %>%
select(-locality) %>%
expand_grid(locality = c(NA_character_, unique(L_dat$locality)))
resp.locality <- newdat.locality %>%
bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>%
pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))
gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) +
geom_point(data = filter(resp.locality.scaled, !is.na(locality)),
mapping = aes(color = locality), cex = 0.5, show.legend = FALSE) +
geom_point(data = filter(resp.locality.scaled, is.na(locality))) +
geom_point(data = filter(L_arith, nin %in% resp.locality.scaled$nin),
color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()
P.abies
Pa_dat <- Filter(function(x)!all(is.na(x)), datlist$P.abies)Pa_arith <- pivot_longer(Pa_dat, 10:ncol(Pa_dat), names_to = "nin", values_to = "are") %>%
group_by(nin) %>%
summarize(daa = sum(are)/10, # units from are to decare
n = sum(wildlings*are)) %>%
mutate(density = n/daa) %>%
arrange(desc(density), desc(n), daa)
scaler <- filter(Pa_arith, nin == contrast) %>%
pull(density)
Pa_arith <- mutate(Pa_arith, fit = density/scaler)| nin | daa | n | density | fit |
|---|---|---|---|---|
| T1 | 13.8 | 34.0 | 2.5 | 2.0 |
| T4 | 689.3 | 843.4 | 1.2 | 1.0 |
| swamp | 33.2 | 28.6 | 0.9 | 0.7 |
| plantation forest | 14.5 | 10.5 | 0.7 | 0.6 |
| T32 | 306.8 | 49.9 | 0.2 | 0.1 |
| fen | 96.8 | 14.6 | 0.2 | 0.1 |
| spring | 6.1 | 0.8 | 0.1 | 0.1 |
| disturbed | 76.0 | 7.9 | 0.1 | 0.1 |
| bog | 11.8 | 1.0 | 0.1 | 0.1 |
| cultivated | 631.3 | 4.9 | 0.0 | 0.0 |
| T2 | 17.7 | 0.1 | 0.0 | 0.0 |
| T6 | 2.9 | 0.0 | 0.0 | 0.0 |
| T31 | 17.7 | 0.0 | 0.0 | 0.0 |
perfectlyseparated <- c('T6', 'T31', 'T2')
removed <- c(perfectlyseparated, excluded)
Pa_dat <- Pa_dat %>%
filter_at(removed, all_vars(. < 0.5)) %>%
filter(!(is.na(seeds.ExP) | is.na(seeds.ExP))) %>%
select(-c('species'))Pa_means <- Pa_dat[, locality.preds] %>% map_dbl(mean)
Pa_sds <- Pa_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
Pa_dat <- modify_at(Pa_dat, .at = i, ~ (. - Pa_means[i])/Pa_sds[i])
}
summary(Pa_dat)## locality wildlings age bio01
## Length:18786 Min. : 0.0000 Min. :-1.4707 Min. :-2.4383
## Class :character 1st Qu.: 0.0000 1st Qu.:-0.7900 1st Qu.:-0.4425
## Mode :character Median : 0.0000 Median :-0.2999 Median : 0.3878
## Mean : 0.0528 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.0000 3rd Qu.: 0.4352 3rd Qu.: 0.6571
## Max. :13.0000 Max. : 1.8238 Max. : 0.8414
## bio19 seeds.ExP seeds.WALD
## Min. :-1.7911 Min. :-0.2017 Min. :-0.579549
## 1st Qu.:-0.3616 1st Qu.:-0.2017 1st Qu.:-0.497256
## Median : 0.1501 Median :-0.2017 Median :-0.371481
## Mean : 0.0000 Mean : 0.0000 Mean : 0.000000
## 3rd Qu.: 1.0679 3rd Qu.:-0.2014 3rd Qu.: 0.002906
## Max. : 1.2141 Max. :14.5311 Max. :10.142177
## relelev plantation forest cultivated disturbed
## Min. :-4.65908 Min. :0.0000000 Min. :0.0000 Min. :0.00000
## 1st Qu.:-0.35233 1st Qu.:0.0000000 1st Qu.:0.0000 1st Qu.:0.00000
## Median : 0.04963 Median :0.0000000 Median :0.0000 Median :0.00000
## Mean : 0.00000 Mean :0.0003784 Mean :0.3357 Mean :0.04047
## 3rd Qu.: 0.51376 3rd Qu.:0.0000000 3rd Qu.:1.0000 3rd Qu.:0.00000
## Max. : 3.39527 Max. :0.4850000 Max. :1.0000 Max. :1.00000
## T32 fen T4 T2
## Min. :0.0000 Min. :0.00000 Min. :0.0000 Min. :0.000000
## 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.000000
## Median :0.0000 Median :0.00000 Median :0.0000 Median :0.000000
## Mean :0.1632 Mean :0.05135 Mean :0.3661 Mean :0.002631
## 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.000000
## Max. :1.0000 Max. :1.00000 Max. :1.0000 Max. :0.499000
## T1 swamp spring bog
## Min. :0.000000 Min. :0.0000 Min. :0.000000 Min. :0.000000
## 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.:0.000000 1st Qu.:0.000000
## Median :0.000000 Median :0.0000 Median :0.000000 Median :0.000000
## Mean :0.003123 Mean :0.0177 Mean :0.003221 Mean :0.005946
## 3rd Qu.:0.000000 3rd Qu.:0.0000 3rd Qu.:0.000000 3rd Qu.:0.000000
## Max. :1.000000 Max. :1.0000 Max. :1.000000 Max. :1.000000
## T6 T31
## Min. :0.0000000 Min. :0.0000000
## 1st Qu.:0.0000000 1st Qu.:0.0000000
## Median :0.0000000 Median :0.0000000
## Mean :0.0005545 Mean :0.0004149
## 3rd Qu.:0.0000000 3rd Qu.:0.0000000
## Max. :0.4980000 Max. :0.4480000paste(names(Pa_dat)[!(names(Pa_dat) %in% c(contrast, removed))], collapse = " + ")## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog"Compare WALD vs ExP
Pa_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog +
(1 | locality)")
Pa_m1 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pa_m1)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 4160.2 4293.5 -2063.1 4126.2 18769
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 1.827 1.352
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.347
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.28264 0.46979 -6.987 2.80e-12 ***
## seeds.WALD 0.44241 0.04472 9.894 < 2e-16 ***
## age -0.42684 0.75175 -0.568 0.570
## bio01 0.23004 0.60744 0.379 0.705
## bio19 0.70574 0.86238 0.818 0.413
## relelev -0.32952 0.06404 -5.145 2.67e-07 ***
## cultivated -6.66521 0.58588 -11.376 < 2e-16 ***
## disturbed -2.88120 0.46898 -6.144 8.07e-10 ***
## T32 -1.39290 0.19489 -7.147 8.86e-13 ***
## fen -1.73769 0.42175 -4.120 3.78e-05 ***
## T1 -0.15064 0.42250 -0.357 0.721
## swamp -1.29798 0.32585 -3.983 6.79e-05 ***
## spring -2.22515 1.54449 -1.441 0.150
## bog -1.09333 1.37612 -0.795 0.427
## seeds.WALD:age 0.05965 0.05143 1.160 0.246
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pa_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev +
cultivated + disturbed + T32 + fen + T1 + swamp + spring + bog +
(1 | locality)")
Pa_m2 <- glmmTMB(Pa_f2, Pa_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pa_m2)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 4275.2 4408.5 -2120.6 4241.2 18769
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 1.861 1.364
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.347
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.08756 0.47254 -6.534 6.41e-11 ***
## seeds.ExP -0.01568 0.06770 -0.232 0.817
## age -0.41037 0.75759 -0.542 0.588
## bio01 0.16883 0.61227 0.276 0.783
## bio19 0.79135 0.86941 0.910 0.363
## relelev -0.36921 0.06050 -6.103 1.04e-09 ***
## cultivated -6.90142 0.60046 -11.494 < 2e-16 ***
## disturbed -2.98567 0.44660 -6.685 2.30e-11 ***
## T32 -1.49012 0.19098 -7.802 6.07e-15 ***
## fen -1.80618 0.41803 -4.321 1.56e-05 ***
## T1 -0.26371 0.41785 -0.631 0.528
## swamp -1.42100 0.32366 -4.390 1.13e-05 ***
## spring -1.77133 1.47632 -1.200 0.230
## bog -1.26941 1.42697 -0.890 0.374
## seeds.ExP:age -0.09854 0.09317 -1.058 0.290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Pa_m1 0 17
## Pa_m2 115 17Continue with WALD.
Compare zero inflation effects
Pa_m3 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ 1)
summary(Pa_m3)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~1
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 4162.2 4303.4 -2063.1 4126.2 18768
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 1.827 1.352
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.347
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.28263 0.46979 -6.987 2.80e-12 ***
## seeds.WALD 0.44241 0.04472 9.894 < 2e-16 ***
## age -0.42684 0.75175 -0.568 0.570
## bio01 0.23004 0.60744 0.379 0.705
## bio19 0.70574 0.86238 0.818 0.413
## relelev -0.32952 0.06404 -5.145 2.67e-07 ***
## cultivated -6.66522 0.58588 -11.376 < 2e-16 ***
## disturbed -2.88121 0.46898 -6.144 8.07e-10 ***
## T32 -1.39290 0.19489 -7.147 8.86e-13 ***
## fen -1.73770 0.42175 -4.120 3.78e-05 ***
## T1 -0.15064 0.42250 -0.357 0.721
## swamp -1.29798 0.32585 -3.983 6.79e-05 ***
## spring -2.22504 1.54444 -1.441 0.150
## bog -1.09325 1.37609 -0.794 0.427
## seeds.WALD:age 0.05965 0.05143 1.160 0.246
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -18.88 2199.10 -0.009 0.993Pa_m4 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ age)
summary(Pa_m4)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~age
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 4082.0 4230.9 -2022.0 4044.0 18767
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 1.645 1.283
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.563
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.70065 0.49367 -3.445 0.000571 ***
## seeds.WALD 0.47071 0.08818 5.338 9.40e-08 ***
## age 1.05118 0.75423 1.394 0.163404
## bio01 0.14233 0.58405 0.244 0.807471
## bio19 0.53330 0.82029 0.650 0.515608
## relelev -0.28260 0.06231 -4.535 5.75e-06 ***
## cultivated -6.66195 0.58627 -11.363 < 2e-16 ***
## disturbed -2.97487 0.47862 -6.216 5.12e-10 ***
## T32 -1.39752 0.19406 -7.201 5.96e-13 ***
## fen -1.87243 0.40895 -4.579 4.68e-06 ***
## T1 -0.12570 0.38750 -0.324 0.745634
## swamp -1.22709 0.30478 -4.026 5.67e-05 ***
## spring -2.36531 1.56487 -1.512 0.130660
## bog -0.99254 1.36192 -0.729 0.466135
## seeds.WALD:age 0.05025 0.10939 0.459 0.645939
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6100 0.2981 2.046 0.0407 *
## age 2.5112 0.2385 10.528 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pa_m5 <- glmmTMB(Pa_f1, Pa_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Pa_m5)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~age + (1 | locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 3787.2 3944.0 -1873.6 3747.2 18766
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2872 0.5359
## Number of obs: 18786, groups: locality, 9
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 5.737 2.395
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.974
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.16094 0.32224 -3.603 0.000315 ***
## seeds.WALD 0.66246 0.09936 6.668 2.60e-11 ***
## age -0.37772 0.47387 -0.797 0.425401
## bio01 -0.40545 0.33732 -1.202 0.229381
## bio19 0.32013 0.44298 0.723 0.469877
## relelev -0.19049 0.07535 -2.528 0.011474 *
## cultivated -6.72680 0.58479 -11.503 < 2e-16 ***
## disturbed -3.40811 0.52840 -6.450 1.12e-10 ***
## T32 -1.61926 0.19840 -8.162 3.31e-16 ***
## fen -1.32823 0.39136 -3.394 0.000689 ***
## T1 -0.10743 0.33060 -0.325 0.745208
## swamp -1.10088 0.27419 -4.015 5.94e-05 ***
## spring -3.90949 1.74967 -2.234 0.025455 *
## bog -1.70301 1.37601 -1.238 0.215850
## seeds.WALD:age 0.01866 0.12540 0.149 0.881723
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6569 0.8732 1.898 0.0578 .
## age 1.0916 0.8401 1.299 0.1938
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Pa_m5 0.0 20
## Pa_m4 294.8 19
## Pa_m1 373.0 17
## Pa_m3 375.0 18m5 is better than m1, m3, and m4: much better AIC, zero-inflation effects
Drop insignificant interaction term?
Pa_f6 <- update.formula(Pa_f1, ~ . - seeds.WALD:age)
Pa_m6 <- glmmTMB(Pa_f6, Pa_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
summary(Pa_m6)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~age + (1 | locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 3785.2 3934.2 -1873.6 3747.2 18767
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.2927 0.5411
## Number of obs: 18786, groups: locality, 9
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 5.711 2.39
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for nbinom2 family (): 0.975
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.15211 0.31936 -3.608 0.000309 ***
## seeds.WALD 0.65391 0.08092 8.081 6.43e-16 ***
## age -0.35620 0.45463 -0.783 0.433335
## bio01 -0.40959 0.33850 -1.210 0.226275
## bio19 0.32207 0.44547 0.723 0.469680
## relelev -0.19114 0.07519 -2.542 0.011021 *
## cultivated -6.73032 0.58456 -11.514 < 2e-16 ***
## disturbed -3.42030 0.52325 -6.537 6.29e-11 ***
## T32 -1.61740 0.19790 -8.173 3.01e-16 ***
## fen -1.32466 0.39052 -3.392 0.000694 ***
## T1 -0.10484 0.33013 -0.318 0.750801
## swamp -1.09958 0.27405 -4.012 6.01e-05 ***
## spring -3.93216 1.74370 -2.255 0.024129 *
## bog -1.71428 1.37360 -1.248 0.212023
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6642 0.8705 1.912 0.0559 .
## age 1.1061 0.8336 1.327 0.1845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Pa_m6 0 19
## Pa_m5 2 20simulationOutput <- simulateResiduals(fittedModel = Pa_m6)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.050344, p-value = 0.016
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0206, p-value = 0.304
## alternative hypothesis: two.sidedDiagnostics look bad — particularly dispersion test.
Pa_m7 <- glmmTMB(Pa_f6, Pa_dat, family = "genpois", ziformula = ~ age + (1 | locality))
summary(Pa_m7)## Family: genpois ( log )
## Formula:
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~age + (1 | locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 3778.7 3927.7 -1870.4 3740.7 18767
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.09272 0.3045
## Number of obs: 18786, groups: locality, 9
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 7.579 2.753
## Number of obs: 18786, groups: locality, 9
##
## Overdispersion parameter for genpois family (): 2.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.86769 0.24637 -7.581 3.44e-14 ***
## seeds.WALD 0.40904 0.04779 8.559 < 2e-16 ***
## age -1.20711 0.37856 -3.189 0.00143 **
## bio01 0.35056 0.23679 1.480 0.13875
## bio19 0.28991 0.36360 0.797 0.42526
## relelev -0.20621 0.07015 -2.940 0.00328 **
## cultivated -6.30372 0.58887 -10.705 < 2e-16 ***
## disturbed -4.04454 0.68253 -5.926 3.11e-09 ***
## T32 -1.47338 0.18570 -7.934 2.12e-15 ***
## fen -1.30504 0.40327 -3.236 0.00121 **
## T1 -0.41010 0.35165 -1.166 0.24352
## swamp -0.98527 0.25745 -3.827 0.00013 ***
## spring -1.95071 1.42452 -1.369 0.17088
## bog -1.12067 1.24805 -0.898 0.36922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.01715 1.18071 -0.014 0.988
## age -0.86046 1.11967 -0.768 0.442simulationOutput <- simulateResiduals(fittedModel = Pa_m7)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.41878, p-value = 0.136
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0141, p-value = 0.464
## alternative hypothesis: two.sidedBetter diagnostics with generalized poisson — including dispersion test.
Pa_m8 <- glmmTMB(Pa_f6, Pa_dat, family = poisson, ziformula = ~ age + (1 | locality))
summary(Pa_m8)## Family: poisson ( log )
## Formula:
## wildlings ~ seeds.WALD + age + bio01 + bio19 + relelev + cultivated +
## disturbed + T32 + fen + T1 + swamp + spring + bog + (1 |
## locality)
## Zero inflation: ~age + (1 | locality)
## Data: Pa_dat
##
## AIC BIC logLik deviance df.resid
## 3908.5 4049.6 -1936.2 3872.5 18768
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 7.368e-10 2.714e-05
## Number of obs: 18786, groups: locality, 9
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 3.031 1.741
## Number of obs: 18786, groups: locality, 9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.008332 0.119825 -0.070 0.944562
## seeds.WALD 0.218690 0.030478 7.175 7.21e-13 ***
## age -0.255051 0.179221 -1.423 0.154706
## bio01 -0.581891 0.126824 -4.588 4.47e-06 ***
## bio19 0.480586 0.160385 2.996 0.002731 **
## relelev 0.126163 0.056431 2.236 0.025372 *
## cultivated -6.693873 0.582473 -11.492 < 2e-16 ***
## disturbed -2.850806 0.435732 -6.543 6.05e-11 ***
## T32 -1.728755 0.168677 -10.249 < 2e-16 ***
## fen -1.068647 0.369595 -2.891 0.003835 **
## T1 -0.149395 0.288148 -0.518 0.604134
## swamp -0.920005 0.241456 -3.810 0.000139 ***
## spring -2.619841 1.395613 -1.877 0.060491 .
## bog -1.659463 1.247526 -1.330 0.183451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0831 0.6029 5.114 3.15e-07 ***
## age 0.8860 0.5953 1.488 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1simulationOutput <- simulateResiduals(fittedModel = Pa_m8)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.65524, p-value = 0.232
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.014, p-value = 0.616
## alternative hypothesis: two.sidedBetter diagnostics with poisson — including dispersion test.
## dAIC df
## Pa_m7 0.0 19
## Pa_m6 6.5 19
## Pa_m8 129.8 18Continue with generalized poisson model: best AIC and OK diagnostics.
Interpret relative establishment probability
selmod <- Pa_m7
terms <- attributes(selmod$frame)$terms
newdat <- attributes(terms)$factors %>%
as_tibble() %>%
mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Pa_means["seeds.WALD"]) / Pa_sds["seeds.WALD"],
age = (Ps_means["age"] - Pa_means["age"]) / Pa_sds["age"],
bio01 = (Ps_means["bio01"] - Pa_means["bio01"]) / Pa_sds["bio01"],
bio19 = (Ps_means["bio19"] - Pa_means["bio19"]) / Pa_sds["bio19"],
relelev = (Ps_means["relelev"] - Pa_means["relelev"]) / Pa_sds["relelev"],
locality = NA_character_) %>%
distinct()
newdat <- newdat %>%
mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>%
replace_na(list(nin = "T4"))Pa_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Pa_link <- mutate(Pa_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>%
select(-se.fit) %>%
mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Pa_link, nin == "T4") %>%
pull(fit)
link.scaled <- mutate_at(Pa_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) +
geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) +
geom_point() +
geom_point(data = filter(Pa_arith, nin %in% link.scaled$nin),
col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()
newdat.locality <- newdat %>%
select(-locality) %>%
expand_grid(locality = c(NA_character_, unique(Pa_dat$locality)))
resp.locality <- newdat.locality %>%
bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>%
pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))
gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) +
geom_point(data = filter(resp.locality.scaled, !is.na(locality)),
mapping = aes(color = locality), cex = 0.5, show.legend = FALSE) +
geom_point(data = filter(resp.locality.scaled, is.na(locality))) +
geom_point(data = filter(Pa_arith, nin %in% resp.locality.scaled$nin),
color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()
P.contorta
Pc_dat <- Filter(function(x)!all(is.na(x)), datlist$P.contorta)Pc_arith <- pivot_longer(Pc_dat, 10:ncol(Pc_dat), names_to = "nin", values_to = "are") %>%
group_by(nin) %>%
summarize(daa = sum(are)/10, # units from are to decare
n = sum(wildlings*are)) %>%
mutate(density = n/daa) %>%
arrange(desc(density), desc(n), daa)
scaler <- filter(Pc_arith, nin == contrast) %>%
pull(density)
Pc_arith <- mutate(Pc_arith, fit = density/scaler)| nin | daa | n | density | fit |
|---|---|---|---|---|
| disturbed | 216.2 | 3603.7 | 16.7 | 25.2 |
| T34 | 117.7 | 233.6 | 2.0 | 3.0 |
| T4 | 578.5 | 383.2 | 0.7 | 1.0 |
| T30 | 2.5 | 1.6 | 0.6 | 1.0 |
| T1 | 27.5 | 11.0 | 0.4 | 0.6 |
| swamp | 4.5 | 1.5 | 0.3 | 0.5 |
| fen | 219.1 | 43.4 | 0.2 | 0.3 |
| plantation forest | 119.1 | 12.1 | 0.1 | 0.2 |
| bog | 90.5 | 3.9 | 0.0 | 0.1 |
| cultivated | 7.6 | 0.0 | 0.0 | 0.0 |
| T27 | 8.0 | 0.0 | 0.0 | 0.0 |
| T32 | 23.3 | 0.0 | 0.0 | 0.0 |
perfectlyseparated <- c('cultivated', 'T27', 'T32')
removed <- c(perfectlyseparated, excluded)
Pc_dat <- Pc_dat %>%
filter_at(removed, all_vars(. < 0.5)) %>%
filter(!(is.na(seeds.ExP) | is.na(seeds.ExP) | is.na(relelev))) %>%
select(-c('species'))Pc_means <- Pc_dat[, locality.preds] %>% map_dbl(mean)
Pc_sds <- Pc_dat[, locality.preds] %>% map_dbl(sd)
for (i in locality.preds) {
Pc_dat <- modify_at(Pc_dat, .at = i, ~ (. - Pc_means[i])/Pc_sds[i])
}
summary(Pc_dat)## locality wildlings age
## Length:12754 Min. : 0.0000 Min. :-1.413778
## Class :character 1st Qu.: 0.0000 1st Qu.:-0.539537
## Mode :character Median : 0.0000 Median : 0.006863
## Mean : 0.3473 Mean : 0.000000
## 3rd Qu.: 0.0000 3rd Qu.: 0.334704
## Max. :131.0000 Max. : 1.755345
## bio01 bio19 seeds.ExP seeds.WALD
## Min. :-1.0009 Min. :-1.0620 Min. :-0.1918 Min. :-0.2757
## 1st Qu.:-0.9670 1st Qu.:-1.0331 1st Qu.:-0.1918 1st Qu.:-0.2751
## Median :-0.8958 Median :-0.6807 Median :-0.1918 Median :-0.2701
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 1.1100 3rd Qu.: 0.8907 3rd Qu.:-0.1918 3rd Qu.:-0.2185
## Max. : 1.1246 Max. : 1.3528 Max. :14.3222 Max. :20.4881
## relelev plantation forest cultivated disturbed
## Min. :-2.2553 Min. :0.000000 Min. :0.000000 Min. :0.0000
## 1st Qu.:-0.4599 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.0000
## Median :-0.1993 Median :0.000000 Median :0.000000 Median :0.0000
## Mean : 0.0000 Mean :0.004167 Mean :0.000744 Mean :0.1665
## 3rd Qu.: 0.2092 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.:0.0000
## Max. : 5.1241 Max. :0.499000 Max. :0.496000 Max. :1.0000
## T32 fen T34 T4
## Min. :0.000000 Min. :0.0000 Min. :0.00000 Min. :0.000
## 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.000
## Median :0.000000 Median :0.0000 Median :0.00000 Median :0.014
## Mean :0.001691 Mean :0.1715 Mean :0.09193 Mean :0.451
## 3rd Qu.:0.000000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:1.000
## Max. :0.495000 Max. :1.0000 Max. :1.00000 Max. :1.000
## T1 swamp T30
## Min. :0.00000 Min. :0.000000 Min. :0.000000
## 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.000000
## Median :0.00000 Median :0.000000 Median :0.000000
## Mean :0.02158 Mean :0.003533 Mean :0.001921
## 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.:0.000000
## Max. :1.00000 Max. :1.000000 Max. :1.000000
## T27 bog
## Min. :0.0000000 Min. :0.0000
## 1st Qu.:0.0000000 1st Qu.:0.0000
## Median :0.0000000 Median :0.0000
## Mean :0.0005656 Mean :0.0707
## 3rd Qu.:0.0000000 3rd Qu.:0.0000
## Max. :0.4940000 Max. :1.0000paste(names(Pc_dat)[!(names(Pc_dat) %in% c(contrast, removed))], collapse = " + ")## [1] "locality + wildlings + age + bio01 + bio19 + seeds.ExP + seeds.WALD + relelev + disturbed + fen + T34 + T1 + swamp + T30 + bog"Compare WALD vs ExP
Pc_f1 <- formula("wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev +
disturbed + fen + T34 + T1 + swamp + T30 + bog + (1 | locality)")
Pc_m1 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pc_m1)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +
## fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Data: Pc_dat
##
## AIC BIC logLik deviance df.resid
## 3424.0 3543.3 -1696.0 3392.0 12738
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.06456 0.2541
## Number of obs: 12754, groups: locality, 6
##
## Overdispersion parameter for nbinom2 family (): 0.0235
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.1381 0.4437 -13.834 < 2e-16 ***
## seeds.WALD 1.1642 0.1759 6.619 3.63e-11 ***
## age 0.1507 0.2331 0.647 0.5178
## bio01 -8.5144 1.2909 -6.596 4.23e-11 ***
## bio19 6.3433 1.3481 4.705 2.53e-06 ***
## relelev -0.1538 0.1986 -0.774 0.4388
## disturbed 3.5411 0.3085 11.479 < 2e-16 ***
## fen 2.5660 1.0236 2.507 0.0122 *
## T34 5.3311 0.9792 5.444 5.20e-08 ***
## T1 6.1689 1.0413 5.924 3.14e-09 ***
## swamp -0.7505 2.5676 -0.292 0.7701
## T30 1.4172 1.7734 0.799 0.4242
## bog -3.4078 0.7741 -4.402 1.07e-05 ***
## seeds.WALD:age -0.5273 0.1211 -4.354 1.34e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pc_f2 <- formula("wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev +
disturbed + fen + T34 + T1 + swamp + T30 + bog + (1 | locality)")
Pc_m2 <- glmmTMB(Pc_f2, Pc_dat, family = "nbinom2", ziformula = ~ 0)
summary(Pc_m2)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.ExP * age + bio01 + bio19 + relelev + disturbed +
## fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Data: Pc_dat
##
## AIC BIC logLik deviance df.resid
## 3458.6 3577.9 -1713.3 3426.6 12738
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.04917 0.2217
## Number of obs: 12754, groups: locality, 6
##
## Overdispersion parameter for nbinom2 family (): 0.022
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.026105 0.424305 -14.202 < 2e-16 ***
## seeds.ExP 0.801158 0.162047 4.944 7.65e-07 ***
## age 0.007757 0.216311 0.036 0.971393
## bio01 -9.267952 1.195267 -7.754 8.91e-15 ***
## bio19 7.309543 1.259363 5.804 6.47e-09 ***
## relelev -0.234139 0.197567 -1.185 0.235972
## disturbed 3.650548 0.310504 11.757 < 2e-16 ***
## fen 2.376049 0.991193 2.397 0.016523 *
## T34 5.458067 0.941805 5.795 6.82e-09 ***
## T1 6.447843 0.998741 6.456 1.08e-10 ***
## swamp 0.008496 2.252186 0.004 0.996990
## T30 1.771512 1.802595 0.983 0.325727
## bog -3.361874 0.775087 -4.337 1.44e-05 ***
## seeds.ExP:age -0.357946 0.104856 -3.414 0.000641 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## dAIC df
## Pc_m1 0.0 16
## Pc_m2 34.6 16Continue with WALD.
Compare zero inflation effects
Pc_m3 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ 1)
summary(Pc_m3)## Family: nbinom2 ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +
## fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Zero inflation: ~1
## Data: Pc_dat
##
## AIC BIC logLik deviance df.resid
## 3426.0 3552.7 -1696.0 3392.0 12737
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.06457 0.2541
## Number of obs: 12754, groups: locality, 6
##
## Overdispersion parameter for nbinom2 family (): 0.0235
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.1381 0.4437 -13.834 < 2e-16 ***
## seeds.WALD 1.1642 0.1759 6.619 3.63e-11 ***
## age 0.1507 0.2331 0.647 0.5179
## bio01 -8.5145 1.2909 -6.596 4.23e-11 ***
## bio19 6.3434 1.3481 4.705 2.53e-06 ***
## relelev -0.1538 0.1986 -0.774 0.4388
## disturbed 3.5411 0.3085 11.479 < 2e-16 ***
## fen 2.5660 1.0236 2.507 0.0122 *
## T34 5.3311 0.9792 5.444 5.20e-08 ***
## T1 6.1689 1.0413 5.924 3.14e-09 ***
## swamp -0.7506 2.5676 -0.292 0.7700
## T30 1.4172 1.7734 0.799 0.4242
## bog -3.4078 0.7741 -4.402 1.07e-05 ***
## seeds.WALD:age -0.5273 0.1211 -4.354 1.34e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.23 2437.67 -0.007 0.995Pc_m4 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ age)
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')Pc_m5 <- glmmTMB(Pc_f1, Pc_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
# 2: In fitTMB(TMBStruc) :
# Model convergence problem; false convergence (8). See vignette('troubleshooting')Overparameterized models? bio01 and bio19 show largest coefficients (5-6) while insignificant for other species. Try dropping them.
Pc_f6 <- update.formula(Pc_f1, ~ . - bio01 - bio19)
Pc_m6 <- glmmTMB(Pc_f6, Pc_dat, family = "nbinom2", ziformula = ~ age + (1 | locality))
# Warning messages:
# 1: In fitTMB(TMBStruc) :
# Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')Pc_m7 <- glmmTMB(Pc_f1, Pc_dat, family = "genpois", ziformula = ~ age + (1 | locality))
summary(Pc_m7)## Family: genpois ( log )
## Formula:
## wildlings ~ seeds.WALD * age + bio01 + bio19 + relelev + disturbed +
## fen + T34 + T1 + swamp + T30 + bog + (1 | locality)
## Zero inflation: ~age + (1 | locality)
## Data: Pc_dat
##
## AIC BIC logLik deviance df.resid
## 3355.0 3496.6 -1658.5 3317.0 12735
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 0.1903 0.4363
## Number of obs: 12754, groups: locality, 6
##
## Zero-inflation model:
## Groups Name Variance Std.Dev.
## locality (Intercept) 6.884e-09 8.297e-05
## Number of obs: 12754, groups: locality, 6
##
## Overdispersion parameter for genpois family (): 87.5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.12923 0.37565 -10.992 < 2e-16 ***
## seeds.WALD 0.30719 0.02676 11.480 < 2e-16 ***
## age 0.36613 0.29156 1.256 0.2092
## bio01 -8.92063 1.57850 -5.651 1.59e-08 ***
## bio19 7.63388 1.53405 4.976 6.48e-07 ***
## relelev -0.17693 0.11906 -1.486 0.1373
## disturbed 3.38649 0.23107 14.655 < 2e-16 ***
## fen 0.71455 0.35603 2.007 0.0447 *
## T34 3.52341 0.37936 9.288 < 2e-16 ***
## T1 5.49293 0.77380 7.099 1.26e-12 ***
## swamp 1.41629 1.62537 0.871 0.3836
## T30 2.52351 1.03456 2.439 0.0147 *
## bog -1.66871 0.65506 -2.547 0.0109 *
## seeds.WALD:age -0.03926 0.02091 -1.878 0.0604 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.161 8.316 -0.500 0.617
## age 2.955 4.737 0.624 0.533Huge overdispersion parameter. Also huge coefficients for bio01 and bio19.
## dAIC df
## Pc_m7 0 19
## Pc_m3 71 17simulationOutput <- simulateResiduals(fittedModel = Pc_m7)## It seems you are diagnosing a glmmTBM model. There are still a few minor limitations associatd with this package. The most important is that glmmTMB doesn't implement an option to create unconditional predictions from the model, which means that predicted values (in res ~ pred) plots include the random effects. With strong random effects, this can sometimes create diagonal patterns from bottom left to top right in the res ~ pred plot. All other tests and plots should work as desired. Please see https://github.com/florianhartig/DHARMa/issues/16 for further details.plot(simulationOutput)
testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted
## vs. simulated
##
## data: simulationOutput
## ratioObsSim = 0.83723, p-value = 0.648
## alternative hypothesis: two.sidedtestZeroInflation(simulationOutput)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.99957, p-value = 0.864
## alternative hypothesis: two.sidedDiagnostics look OK.
Continue with generalized poisson model: best AIC and OK diagnostics.
Interpret relative establishment probability
selmod <- Pc_m7
terms <- attributes(selmod$frame)$terms
newdat <- attributes(terms)$factors %>%
as_tibble() %>%
select(-`seeds.WALD:age`) %>%
mutate(seeds.WALD = (Ps_means["seeds.WALD"] - Pc_means["seeds.WALD"]) / Pc_sds["seeds.WALD"],
age = (Ps_means["age"] - Pc_means["age"]) / Pc_sds["age"],
bio01 = (Ps_means["bio01"] - Pc_means["bio01"]) / Pc_sds["bio01"],
bio19 = (Ps_means["bio19"] - Pc_means["bio19"]) / Pc_sds["bio19"],
relelev = (Ps_means["relelev"] - Pc_means["relelev"]) / Pc_sds["relelev"],
locality = NA_character_) %>%
distinct()
newdat <- newdat %>%
mutate(nin = names(newdat)[apply(newdat, 1, match, x = 1)]) %>%
replace_na(list(nin = "T4"))Pc_link <- bind_cols(newdat, predict(selmod, newdat, se.fit = TRUE, type = "link"))
Pc_link <- mutate(Pc_link, ci.l = fit - 1.96*se.fit, ci.u = fit + 1.96*se.fit) %>%
select(-se.fit) %>%
mutate_at(c('fit', 'ci.l', 'ci.u'), ~ exp(.))
scaler <- filter(Pc_link, nin == "T4") %>%
pull(fit)
link.scaled <- mutate_at(Pc_link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))gg <- ggplot(link.scaled, aes(reorder(nin, fit), fit)) +
geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u)) +
geom_point() +
geom_point(data = filter(Pc_arith, nin %in% link.scaled$nin),
col = "grey", pch = 1)
gg + scale_y_log10() + coord_flip()
newdat.locality <- newdat %>%
select(-locality) %>%
expand_grid(locality = c(NA_character_, unique(Pc_dat$locality)))
resp.locality <- newdat.locality %>%
bind_cols(fit = predict(selmod, newdat.locality, type = "response"))
scaler <- filter(resp.locality, nin == "T4" & is.na(locality)) %>%
pull(fit)
resp.locality.scaled <- mutate_at(resp.locality, 'fit', ~ `/`(., scaler))
gg <- ggplot(resp.locality.scaled, aes(reorder(nin, fit), fit)) +
geom_point(data = filter(resp.locality.scaled, !is.na(locality)),
mapping = aes(color = locality), cex = 0.5, show.legend = FALSE) +
geom_point(data = filter(resp.locality.scaled, is.na(locality))) +
geom_point(data = filter(L_arith, nin %in% resp.locality.scaled$nin),
color = "grey50",pch = 1)
gg + scale_y_log10() + coord_flip()## Warning: Transformation introduced infinite values in continuous y-axis
All species
link.list <- list(Ps_link, L_link, Pa_link, Pc_link)
names(link.list) <- c('Picea sitchensis/lutzii', 'Larix', 'Picea abies', 'Pinus contorta')
link <- bind_rows(link.list, .id = 'species')
scaler <- filter(link, species == 'Picea sitchensis/lutzii', nin == "T4") %>%
pull(fit)
link.scaled <- mutate_at(link, c('fit', 'ci.l', 'ci.u'), ~ `/`(., scaler))arith.list <- list(filter(Ps_arith, nin %in% Ps_link$nin),
filter(L_arith, nin %in% L_link$nin),
filter(Pa_arith, nin %in% Pa_link$nin),
filter(Pc_arith, nin %in% Pc_link$nin))
names(arith.list) <- c('Picea sitchensis/lutzii', 'Larix', 'Picea abies', 'Pinus contorta')
arith <- bind_rows(arith.list, .id = 'species')
scaler <- filter(arith, species == 'Picea sitchensis/lutzii', nin == "T4") %>%
pull(density)
arith.scaled <- mutate(arith, fit = density/scaler)link.scaled[,c('species','nin','fit','ci.l','ci.u')]## # A tibble: 45 x 5
## species nin fit ci.l ci.u
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 Picea sitchensis/lutzii T4 1 0.599 1.67
## 2 Picea sitchensis/lutzii cultivated 0.0719 0.0432 0.120
## 3 Picea sitchensis/lutzii disturbed 0.466 0.268 0.811
## 4 Picea sitchensis/lutzii T32 0.448 0.264 0.760
## 5 Picea sitchensis/lutzii fen 0.853 0.510 1.43
## 6 Picea sitchensis/lutzii T34 1.03 0.618 1.72
## 7 Picea sitchensis/lutzii T13 1.60 0.250 10.3
## 8 Picea sitchensis/lutzii T2 2.85 1.42 5.70
## 9 Picea sitchensis/lutzii T1 0.0974 0.0502 0.189
## 10 Picea sitchensis/lutzii swamp 0.131 0.0625 0.276
## # … with 35 more rowsgg <- ggplot(link.scaled, aes(fct_reorder(nin, fit, first), fit, color = species)) +
geom_linerange(mapping = aes(ymin = ci.l, ymax = ci.u),
position = position_dodge(width = 0.5)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_point(data = arith.scaled, pch = 1, position = position_dodge(width = 0.5))
gg + scale_y_log10() + coord_flip() +
labs(x = 'NiN type', y = 'relative establishment probability', title = 'Relative establishment across species and NiN types')
Filled points are model estimates with lines showing 95% confidence intervals. Open points are calculated directly from the number of wildlings per unit area.
sessionInfo
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] DHARMa_0.2.4 glmmTMB_0.2.3 forcats_0.4.0 stringr_1.4.0
## [5] dplyr_0.8.3 purrr_0.3.2 readr_1.3.1 tidyr_1.0.0
## [9] tibble_2.1.1 ggplot2_3.1.1 tidyverse_1.2.1 here_0.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.0 jsonlite_1.6 splines_3.6.0
## [4] foreach_1.4.4 modelr_0.1.4 gap_1.2.1
## [7] assertthat_0.2.1 highr_0.8 stats4_3.6.0
## [10] cellranger_1.1.0 yaml_2.2.0 numDeriv_2016.8-1.1
## [13] pillar_1.4.0 backports_1.1.4 lattice_0.20-38
## [16] glue_1.3.1 bbmle_1.0.22 digest_0.6.19
## [19] rvest_0.3.4 minqa_1.2.4 colorspace_1.4-1
## [22] htmltools_0.4.0 Matrix_1.2-17 plyr_1.8.4
## [25] spaMM_3.0.0 pkgconfig_2.0.2 broom_0.5.2
## [28] haven_2.1.0 mvtnorm_1.0-10 scales_1.0.0
## [31] lme4_1.1-21 proxy_0.4-23 generics_0.0.2
## [34] ellipsis_0.3.0 withr_2.1.2 pbapply_1.4-2
## [37] TMB_1.7.15 lazyeval_0.2.2 cli_1.1.0
## [40] magrittr_1.5 crayon_1.3.4 readxl_1.3.1
## [43] evaluate_0.13 fansi_0.4.0 nlme_3.1-139
## [46] MASS_7.3-51.4 xml2_1.2.0 tools_3.6.0
## [49] hms_0.4.2 lifecycle_0.1.0 munsell_0.5.0
## [52] compiler_3.6.0 rlang_0.4.2 grid_3.6.0
## [55] nloptr_1.2.1 iterators_1.0.10 rstudioapi_0.10
## [58] rmarkdown_1.12 boot_1.3-22 gtable_0.3.0
## [61] codetools_0.2-16 R6_2.4.0 lubridate_1.7.4
## [64] knitr_1.23 bdsmatrix_1.3-4 zeallot_0.1.0
## [67] utf8_1.1.4 rprojroot_1.3-2 stringi_1.4.3
## [70] parallel_3.6.0 Rcpp_1.0.1 vctrs_0.2.1
## [73] tidyselect_0.2.5 xfun_0.7